Superconducting states and depinning transitions of Josephson ladders 
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We present analytical and numerical studies of pinned su- 
perconducting states of open-ended Josephson ladder arrays, 
neglecting inductances but taking edge effects into account. 
Treating the edge effects perturbatively, we find analytical ap- 
proximations for three of these superconducting states - the 
no- vortex, fully-frustrated and single- vortex states - as func- 
tions of the dc bias current I and the frustration /. Bifurca- 
tion theory is used to derive formulas for the depinning cur- 
rents and critical frustrations at which the superconducting 
states disappear or lose dynamical stability as I and / are var- 
ied. These results are combined to yield a zero-temperature 
stability diagram of the system with respect to I and /. To 
highlight the effects of the edges, we compare this dynami- 
cal stability diagram to the thermodynamic phase diagram 
for the infinite system where edges have been neglected. We 
briefly indicate how to extend our methods to include self- 
inductances. 

PACS Numbers: 74.50. +r, 05.45.+b, 74.40. +k 



I. INTRODUCTION 

Arrays of Josephson Junctions are of interest in sev- 
eral branches of physical . They have many technological 
applications, including high-frequency emitters and de- 
tectors, parametric amplifiers, local oscillators, and volU 
age standardsEra. They also shed light on the structurally 
and pinningu properties of the high-T c superconducting 
cuprates. At the same time, they provide model sys- 
tems for the study of -problems in both spatiotempo- 
ral nonlinear, dynamicaTB and nonequilibrium statisti- 
cal physicsEM. For instance, the depinning transitions 
and nonlinear wave propagation seen in Josephson ar- 
rays are analogous to those found in incommensurate 
systems, earthquake models, Type II superconductors, 
and charge-density waves. 
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From the standpoint of dynamical systems theory0, 
Josephson arrays can be viewed as large collections of 
coupled nonlinear oscillators. Unfortunately, because 
of their nonlinearity and large number of degrees of 
freedom, these arrays are inherently difficult to analyze 
mathematically. A further complication is that there is 
an intrinsic physical coupling among junctions, due to 
fluxoid quantization, which is more awkward to handle 
than the nearest-neighbor interaction usually assumed in 
idealized models of coupled oscillators. And when the 
effects of self-fields and inductances are included, there 
is even less hope of making analytical progress. 

Despite these obstacles, some encouraging advances 
have occurred recently in the mathematical analysis of 
Josephson arrays, especially for one-dimensional (1 
systems where the junctions are connected in serie; 
or in parallel! The logical next step is to tackle two- 
dimensional (2-D) arrays. Much of the previous theoret- 
ical work on 2-D arrays has focused on numerical simula- 
tion of the current-voltage characteristicaij, in an effort 
to link the rich spatiotemporal dynamics of 2-D arrays to 
the averaged quantities that are most readily measured 
experimentally. On the mathematical side, there are re- 
cent indications that 2-D arrays, like their 1-D counter- 
parts, are also going to be tractable in some regimesE3. 

An ideal example to explore the crossover between 1- 
D and 2-D behavior is tbe._k>sephson ladder array (Fig- 
ure Q). Following Kardanijii-J, several authors have stud- 
ied various statistical properties of the frustrated ladder, 
including its ground state, the complicated landscape of 
solutions at zero temperaturep-the, low- lying excitations, 
and the linear response regimalTO. However, all of these 
authors restricted attention to ladders in the absence of 
a driving current. Only recently has the fully dynamical 
problem been addressed, through numerical simulatkuas 
of the depinning transitionE£l and vortex propagationEM 

In this paper we use the tools of nonlinear dynam- 
ics to analyze the superconducting states of ladder ar- 
rays, p^lther dynamical regimes will be discussed else- 
whereEJEil.) Mathematically, the superconducting states 
correspond to fixed points of the governing circuit equa- 
tions. We study the bifurcations of these fixed points 
with respect to variations in the applied dc bias current 
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I and the frustration / introduced by an external mag- 
netic field. As / is increased from at fixed /, we find 
that the stable superconducting states are destroyed in 
saddle-node bifurcations at certain critical values of the 
current. Then the system depins from its original static 
configuration and evolves toward some other state. The 
new state might be another fixed point, or it might be 
a running solution, in which case a nonzero dc voltage 
appears across the array. Global depinning of the ar- 
ray occurs when the last stable fixed point is destroyed. 
One of the main results of our analysis is a set of ap- 
proximate analytical formulas for the critical currents at 
which the fixed points are destroyed, as a function of /, 
for the three most important types of superconducting 
states: the no-vortex, fully-frustrated and single-vortex 
solutions. 

Another important finding is that symmetry plays a 
crucial role in the dynamics of the ladder. As we will 
show below, much of the behavior of the ladder can be 
understood by focusing on states that are "up-down sym- 
metric" - in other words, states where the phases of the 
horizontal junctions on the top and bottom of any given 
plaquette are equal in magnitude but opposite in sign 
at all times. All of the stable superconducting states 
possess this symmetry. Even when the parameters are 
chosen so that depinning occurs, the subsequent tran- 
sients and long-term running solutions typically remain 
up-down symmetric. But there is at least one excep- 
tion: when the single-vortex state is destabilized by low- 
ering the frustration / below some critical value / m in, 
the system depins via a symmetry-breaking bifurcation. 
During the transient behavior, the up-down symmetry is 
lost temporarily, but is then recovered as the system ex- 
pels flux from the array and evolves toward the no- vortex 
state. 

These symmetry considerations establish an unex- 
pected link between the study of Josephson arrays and 
some recent developments in nonlinear dynamics. In 
mathematical terms, the up-down symmetric states of a 
ladder array form an invariant manifold of the full state 
space. A symmetry-breaking bifurcation occurs when 
this manifold loses stability in a transverse direction. The 
same issue - the transverse stability of an invariant man- 
ifold - arises in the study of riddled basins, synchronized 
chaos, on-off intermittency, and blowout bifurcationscj. 
These connections suggest a promising line of future re- 
search on Josephson arrays, particularly with regard to 
their chaotic states. 

We will also show that the dynamics of ladder arrays 
are very strongly influenced by edge effects. One might 
have supposed these effects to be negligible, especially 
in long ladders, since their influence on the supercon- 
ducting solutions dies off exponentially fast away from 
the boundaries. Yet although the edges do indeed have a 
small effect on the form of the superconducting solutions, 
they have a large effect on the stability of those solutions. 
Much of this paper is devoted to investigating the effects 
of the edges, first on the superconducting states them- 



selves, then on their stability, and finally on the entire 
phase diagram. 

This paper is organized as follows. Section |l| reviews 
the model equations for the ladder and discusses their 
symmetry properties. In Section III we obtain analyti- 
cal approximations for three numerically observed super- 
conducting solutions: the no- vortex (nv), single- vortex 
(sv), and fully-frustrated (ff) solutions. In all cases, edge 
effects are taken into account via perturbation theory. 
Next, in Section |iy| we describe the dynamical simula- 
tions which reveal the depinning properties of the nv, 
sv, and ff configurations and relate them to the global 
depinning of the ladder. In Section [V| we establish the 
rigorous connection of the dynamical depinning with the 
stability of these three fixed points. For all of them, 
we characterize the bifurcations and study their stabil- 
ity diagrams, both in the presence and in the absence of 
edges. When possible, analytical approximations to the 
critical currents are obtained. We show that the depin- 
ning transitions correspond to saddle-node bifurcations 
that are edge-dominated for almost all values of the frus- 
tration. Moreover, we find that some of the supercon- 
ducting states can be destabilized via a subcritical pitch- 
fork bifurcation as the frustration is reduced; in physical 
terms, this is a symmetry-breaking bifurcation in which 
flux is expelled transver sally from the ladder. In the fi- 
nal section we summarize our conclusions in two phase 
diagrams (for ladders with and without edges), and we 
relate our results to those found by previous authors. 
We also add two more technical appendices: Appendix 
A compares the single-vortex configuration in the ladder 
with the corresponding kink-like solution in 1-D parallel 
arrays; Appendix B briefly indicates how to extend our 
approach to include self- inductance effects. 



II. THE SYSTEM 

We study an open-ended Josephson ladder with N 
square plaquettes, i.e., an array formed by two rows of 
N+l weakly coupled superconducting islands (Figure [j]). 
The array is driven by a perpendicular uniform dc current 
/, and a magnetic field is applied transverse to the plane 
of the device. Each weak link between islands constitutes 
a junction. Its state is described by the gauge-invariant 
phase difference (f>j , arising from the macroscopic charac- 
ter of the quantum wavefunction of the superconductors. 

Assuming zero temperature, negligible charging (quan- 
tum) effects, and identical junctions, the dynamics af 
each junction is given, in the three-channel RCSJ modelEJ, 
by the nonlinear differential equation 



(1) 



Here (3 C is the McCumber parameterEJ; Ij is given m 
units of the critical current of each (identical) junction; 
derivatives are with respect to time normalized in units 
of ijJj 1 , the inverse of the plasma frequency; and N is 
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In summary, the governing equations can be written 
compactly as f(x) = 0, with 

x = (^,...,0^ +1 ,0f,...,^,^,...,^), and with 
the components of f(x) defined by 

£,-(x) = I + M(<j>f_ 1 )-M(<j>f) -Af(4>J), 

j = l,...,N + l (12) 
{ N+1+j ( X )=Af(<t>f)+Af(<j>>?), 

j = l,...,N (13) 

tw+i+i(x) = <t>] + $ - 4>J+x - 4>f - Mnj - f), 

j = l,...,N. (14) 



The dynamical evolution of the state {x(t),x(i)} of the 
system is obtained by numerically solving this system of 
coupled differential and algebraic equations. The dynam- 
ics also depend implicitly on the parameters /, /3 C , /, and 
N. 

In experiments, the most convenient way to probe 
the dynamics of the array is to measure its dc current- 
voltage (IV) characteristics. From the Josephson rela- 
tionsa, the time-dependent voltage across each junction 
is directly proportional to <f>j, the time-derivative of its 
phase. Hence, the total dc voltage V across the array in 
the vertical direction is proportional to the spatial and 
temporal average of all the phase derivatives. Although 
because of this averaging a great deal of dynamical in- 
formation is lost about the spatiotemporal state of the 
system, the IV curve still provides a useful (if somewhat 
coarse) indicator of changes in the underlying dynamics 
as the drive current / is varied. 

In the case of ladder arrays, the IV curves display 
three,pjgions associated with distinct dynamical behav- 



iorsE3E3. At low /, the system is superconducting (V — 0) 
with pinned or slightly oscillating junctions (in the me- 
chanical analog, the pendula are at rest or librating 
slightly). At a depinning current Idep, the array jumps to 
the flux-flow region, in which a finite voltage is produced 
by vortices of magnetic flux moving across the array. At 
still higher currents, the dynamics is characterized by 
whirling modeaj (in the mechanical analog, all the pen- 
dula rotate over the top at a nearly uniform angular ve- 
locity proportional to the applied torque). In this state, 
there is a linear ohmic dependence of V on /. In the 
remainder of the article we focus on the superconducting 
states and the critical current Jd cp at which depinning 



III. OBSERVED SUPERCONDUCTING 
SOLUTIONS AND ANALYTICAL 
APPROXIMATIONS 

We have performed dynamical simulations of the array 
in which the current is ramped up from zero with differ- 
ent initial conditions. The superconducting solutions ob- 
served in those simulations are always static states, i.e., 
fixed points of the system. (In principle, time-dependent 



solutions with high cancelling symmetry could also have 
zero total dc voltage, but we never see such states in 
our simulations.) More specifically, for any given (3 C and 
N , only three types of configurations appear in the nu- 
merics: no-vortex solutions (Figure ^) for / smaller than 
~ 0.3, and solutions of the single-vortex type (Figure |^) 
and fully-frustrated type (Figure ||) for / — > 1/2. Far 
from the edges, the no-vortex (nv) state is characterized 
by identical phases for all junctions. The same applies 
to the single-vortex (sv) configuration far from both the 
edges and the center of the vortex. On the other hand, 
the fully- frustrated (ff) state has a spatial oscillatory pat- 
tern with a wavelength equal to two plaquettes. All 
of these states are modified by noticeable edge effects. 
Although there are many other static solutions of the 
system, our numerical simulations indicate that the no- 
vortex, single-vortex, and fully-frustrated states are the 
only ones needed to explain the depinning behavior of 
the array. 

When dealing with fixed points with up-down symme- 
try ( |Io| ) the defining equations (|l2"|)~(14) become 



kH 



J 



<f>Y +1 -2<ftf=2n(n j -f), j 



,N + 1 (15) 
,N (16) 



where we have defined artificial phases <j>Q — 4>n+i = 0- 
Note that the McCumber parameter (3 C does not appear 
in the equations for the fixed points, and hence does not 
affect their existence. This is consistent with the numer- 
ically observed independence of the depinning behavior 
on C . 

In the rest of this section we obtain analytical approx- 
imations for the nv, sv, and ff configurations mentioned 
above. We follow a common scheme for all of them. 
First, we obtain a no-edge approximation (denoted with 
a dagger f ) for the infinite ladder. Second, we introduce 
the effect of the edges perturbatively to obtain an edge- 
corrected approximation (identified by a double dagger 
X). The calculated configurations have been exhaustively 
compared with the results of numerical simulations with 
excellent agreement. 



A. No-vortex solution 

Figure || (a)-(b) shows a plot of the no-vortex solu- 
tion, as computed numerically, along with the analytical 
approximation described below. This state is character- 
ized by the absence of topological vortices (rij — TVj+i = 

0, Vj) and, far from the edges, by the constancy of the 
phases. 

As a first approximation, let {<fiY ,(f>j* } denote the 
phases of the no- vortex solution for the infinite ladder, 

1. e., in the absence of edge effects. To ease the notation, 
let <f) H ^ denote the common phase of the horizontal junc- 



tions (so 



for all j), and define tp vj< similarly 
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for the phases of the vertical junctions. The only physi- 
cally acceptable solution of equations (|l5])-([l6|) that also 
satisfies the stability condition (Jll|) is: 



arcsin/, 



,ift 



7T/, 



(17) 



where < / < 1/2 and all the angles are restricted to 
the first quadrant. 

This solution exists if and only if I < 1 . A linear stabil- 
ity analysis shows that, for all I < 1, the solution is stable 



if 



= arcsin/ 6 [0,7r/2). The other possible solutions 

with <f> v l—^= it — arcsin/ or (f> H ^ — nf — n are linearly 
unstableEl. In summary when the edges are completely 
neglected, the array behaves like a single junction: its 
only stable no- vortex solution of the observed form ( |l7| ) 
disappears at I — 1 through a saddle-node bifurcation. 
This existence criterion will be used in Section [v] when 
discussing the stability properties of the nv solution. 

Figure |2j shows that this infinite-ladder approximation 
works well near the center of the ladder, but breaks down 
close to the edges. We now take edge effects into account 
by considering an edge-corrected solution (denoted by J) 



(18) 



where {Aj, Bj} denote the corrections. From the fixed 
point equations ([l5|), jl^), the {Aj,Bj} must satisfy 

J + sin(7r/ - = 

sin(7r/ — Bj) + sin(arcsin/ + Aj) (19) 



Aj - A j+1 + 2Bj 



0. 



(20) 



The corrections {Aj , Bj} are expected to be small, except 
in a region very close to the edges. Thus, Equation jl9| ) 
can be expanded to first order in Aj and Bj to obtain a 
second-order difference equation for Aj : 



Aj+i - 2aAj + Aj-! = 

with a = 1 H — , 

cos it J 

from which the horizontal corrections are 

_ A j+1 - Aj 
3 2 

The general solution of (^l|) is 

Aj = p r ^-^ N+1 '> +Qr l - j 
where r — a + \/ a 2 — 1 = e 1 ^. 



(21) 
(22) 

(23) 



(24) 
(25) 



Hence, the edges produce corrections that decay expo- 
nentially from both ends with a characteristic length 
X(I,f). This A = 1/lnr is a measure of how small per- 
turbations decay inside of a region with the no-vortex 
superconducting solution. A similar result was recently 
obtained by Denniston and TangcH using the transfer ma- 
trix method for the particular 1 = case. 
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FIG. 2. No-vortex superconducting solution, (a) Phases 
of the horizontal junctions for a 25 x 1 array with I = 0.25 
and / = 0.2: numerically observed solution (+), and approx- 
imate solution <pf (o). (b) Same as in (a) for the vertical 
junctions. The predicted phases of the infinite-ladder solu- 
tion {(f)j\(f>f^} (short-dashed line) are seen to be valid far 
from the edges in both (a) and (b). (c) Exponential decay 
of the correction from the edges with varying magnetic field 
/ and current I. The symbols correspond to the numerically 
calculated vertical phases. The solid lines come from our ap- 
proximate solution <j}J^ and have slopes ±lnr given by 

To complete the solution, the constants P, Q in fl24| ) 
have to be fixed from the boundary conditions 

I = sin(7r/ - B x ) + sin(arcsin I + A x ) (26) 
I + sm(nf — Bn) = sin(arcsin/ + Ajf+i), (27) 

which result from current conservation at nodes 1 and 
iV + 1 respectively. Since {Aj, Bj} become largest at the 
edges, equations (|2^)-(|27|) are solved numerically without 
linearization. When the array is long enough, such that 
A <C iV+l, the effect of one edge on the other is negligible 
and the solution is further simplified as equations (|2q)- 
( p7| ) decouple. Then Q and P are obtained independently 
by solving 



I = sin ( nf 



Q 



+ sin(arcsin/ + (3) (28) 



5 



I + sin I 7r/ 



P 
~2 



1 - 



= sin(arcsin/ + P). (29) 



Figure g (a)-(b) shows that the above analytical solu- 
tion agrees well with the results of simulations for long 
(N = 25) arrays. The approximation accounts well for 
the effect of the open ends also for short (N = 7) ar- 
rays (not shown). The exponential decay of the per- 
turbation from the edges is also checked satisfactorily 
in Figure ||(c). As expected on physical grounds, the 
edge effects become more important as both the field 
and the current are increased. Thus, the approximation 
is best when / and / are small, and worsens as I — > 1 and 
/ — >l/2(/ = l/2isa singular limit, as seen from the 
vanishing denominator in (|22|)). This establishes limits 
on the use of this approximation for the prediction of the 
depinning current at high values of the frustration. 

Our analytical approximation also explains other fea- 
tures of the numerics. For instance, the corrections 
Aj , Bj are spatially asymmetric with respect to the cen- 
ter of the array when I > — as seen in Figure ||(a)-(b) 
by comparing the rightmost and leftmost phases. Note 
also that for < / < 1/2 the largest vertical phase oc- 
curs at the right end of the ladder — in obvious connec- 
tion with the preferred direction for flux-propagation in 
the array (— x). Moreover, it can readily be shown that 
the change of the frustration from / to 1 — / has only one 
effect: the vertical phases for frustration 1 — / are a mir- 
ror image, with respect to the center of the array, of the 
vertical phases with frustration /. This implies that the 
depinning current will be identical for both values of the 
frustration—as expected, but the direction of propagation 
is reversedEl 



B. Single- vortex solution 

An analytical approximation for the single- vortex con- 
figuration can be obtained in a similar fashion by real- 
izing that the effect of a vortex located in cell a of the 
array is similar to the edge effects in the no- vortex state. 
Note how, if the phases in Figure ||(a)-(b) were reduced 
to [— 7T, 7r), the single- vortex configuration is composed 
of two halves, each of which is equivalent to a no- vortex 
superconducting solution when we move away from the 
edges and from the vortex center a. 

Hence, the zeroth order single-vortex solution 
\4>Y ,<t>? } is identical to the no-edge nv solution given 

w 



m yiY|j. And the edge and vortex-corrected approxima- 
tion with a vortex distribution n a 
a is given by 



1 and rij = 0, Vj 7^ 



Ht 



arcsin I + Aj, <jA = nf — Bj , 



(30) 



where the corrections {Aj, Bj} result now both from the 
presence of the edges and of the vortex in plaquette a. 
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FIG. 3. Single- vortex solution, (a) Phases of the horizontal 
junctions for a 25 x 1 ladder for I — 0.1 and / = 0.2 with a 
vortex in the central plaquette: numerical configuration from 
dynamical simulations (+), and analytical approximation (o) 
as given by (^TJ). (b) Same as in (a) for the vertical junc- 
tions. The phases are not reduced to the [— vr, n) interval, (c) 
Exponential decay of the correction Aj from both the edges 
and from the center of the vortex for different / and I. Sym- 
bols represent numerical simulations and the solid lines are 
the corresponding analytical predictions given by (^). The 
different graphs have been offset for clarity. 

Following the same steps as for the nv configuration, we 
obtain identical expressions (|l9|)-(|25|) for the corrections 
for each half of the array: 



An 



Pr j - a + Qr 1 ~ j , j<a 

p, r j-{N+l) + Q' r a+l-j : j >a 
Aj+1 - Aj 



(31) 
(32) 



where r is, once more, given by (p5|). Therefore, the 
single- vortex solution is obtained by matching two edge- 
corrected no- vortex solutions. In fact, the vortex in cell 
a effectively introduces two new "edges" , at a and a + 1 , 
which also produce similar exponentially decaying correc- 
tions. The matching condition at a and a + 1 is given by 
the fluxoid quantization condition in the cell containing 
the topological vortex: A a — A a+i + 2B a = 2ir. Thus, 



G 



B a 



p' r a 



-N 



Q'-P 



(33) 



This completes the equations needed to determine the 
unknowns P,Q,P',Q' i-u. our solution (|3l|). They can be 
calculated numeric allyE3, for given / and /, from j3^ ) 
and the boundary conditions from current conservation 
at nodes 1, a, a + 1 and N + 1. The approximation is 
compared with numerical simulations in Figure |^ with 
excellent agreement, especially at small /. 

We also note that although the ladder equations (|i~5|)- 
( |l6| ) can be reduced, .-approximately to a discrete sine- 
Gordon equationli3li-lt2l, our analytical expression ([SO]) is 
a better approximation than the much-used sine-Gordon 
kink, which is a good description in strictly 1-D parallel 
arrayscS. A detailed comparison of both approaches is 
presented in Appendix |A|. 



C. Fully- frustrated solution 

The other relevant superconducting state is the fully- 
frustrated solution, which appears in simulations when 
/ w 1/2 (Figure |^). To obtain an analytical approxi- 
mation, we follow once more the same scheme as above: 
first, calculate a no-edge basic solution; then, introduce 
the edges perturbatively. 

In this case, the basic solution is seen numerically to 
oscillate in space with a wavelength equal to two plaque- 
ttes. Thus, when edges are neglected (the infinite-ladder 
approximation) , the phases can be approximated in gen- 
eral by 



4>y* = 2tt [o+(-iy&] 
0f f = 27r [c+(-iyd] 



(34) 



where a, 6, c, d are to be determined from (|l5|)-- ([l6| ) with 
fij = [1 T (— 1) 3 ]/2- First, substitution in Eq. ( [hq ) gives 
c and d 

c = //2-l/4 
d = 6±l/4. 

Second, from Eq. (|l6|) we obtain 



(35) 
(36) 



sin 2na cos 2-Kb = I 
sin 2Kb cos 2ira — 2 sin tt/ cos 27r6 = 0, 

from which we then solve explicitly for a and b in terms 
of the parameters / and / 



1 

27 



b = — arccos J21 2 IL 
2tt v 1 



where 



L = (1 + I 2 ) ± J (1 - I 2 ) 2 - 16I 2 sin 2 nf. 



(37) 
(38) 

(39) 



Figure || (a)-(b) compares the analytical infinite-ladder 
fully-frustrated approximation (|34j)-(j39"|) with numerical 
simulations. The agreement is good except near the ends, 
as expected. 
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FIG. 4. Fully-frustrated solution, (a) Phases of the hori- 
zontal junctions for a 25 x 1 with / = 0.23 and / = 0.5. The 
numerical solution (+) and infinite-ladder analytical approx- 
imation fi- ( — o — ) are seen to coincide except close to the 
edges. The solid line (a guide to the eye) emphasizes the wave- 
length equal to two plaquettes. (b) Same as (a) for the verti- 
cal junctions, (c) Exponential decay of the corrections from 
the edges in the fully-frustrated solution for varying / and I. 
Both the odd (C" um ) and even (Df um ) site corrections have a 
characteristic length ±(lnr)/2 given by (fis|). The solid lines 
are best linear fits with slopes —0.805 (for / = 0.4, / = 0.2) 
and 0.502 (for / = 0.5, / = 0.23). They are in excellent agree- 
ment with the predicted values (lnr)/2 from ([til), which are 
0.806 and 0.488, respectively. 

The approximation above also yields an existence crite- 
rion for the no-edge fully- frustrated solution. In (|39"|), the 
expression inside the square root must be non-negative; 
hence the infinite-ladder fully-frustrated solution does 
not exist if 



I>Ib. 



th 



4 sin nf + 1 - 2 sin nf, 



(40) 



where the subscript "ff" denotes fully-frustrated and "th" 
denotes a theoretical approximation of a bifurcation con- 
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dition. We will use this condition (|40|) later when we 
discuss the depinning of the fully-frustrated solution. 

This predicted form of the infinite-ladder fully- 
frustrated solution agrees with previous findings obtained 
for the special case when there is no driving currentlia. 
For 1 = 0, out solution (|34|)-(fs9l) reduces to the stable 
configurationEa 



sin 2 nf + cb 2 cos 2nf - 2 [(sin 2 nf - sa 2 )sb 2 + ca 2 ] 



(fif 1 = (-l)J' arctan(2), <j>f f = arctan(l/2) 



which coincides with the groipd state for / = 1/2 and 
1 = calculated by BenedicttLj. (To obtain this result 
from the expressions above, note that 2I 2 /L — > (1 + 
4sin 2 7r/) _1 as I — > 0, for the solution corresponding to 
the minus sign in (|39|).) We also note that the infinite- 
ladder fully-frustrated solution exists for all f when I = 
0, i.e., there is no critical magnetic field below which it 
ceases to exist, although it is energetically most favorable 
when / fa 1/2. 

The physical meaning of this solution is clear: it con- 
tains a topological vortex in every other cell, as seen from 
the alternating sequence of zeros and ones for the plaque- 
tte integers {rij} in the equations (p|)-(|l6|). 

In fact, although the solutions with {n oc id = lj n ovcn — 
0} and {n dd = 0, n cvcn = 1} are degenerate in an in- 
finite array, they are not so if the array is finite. How- 
ever, we will show in Section [v| that the depinning of the 
fully-frustrated state is basically unaffected by the parity 
of the number of cells in the ladder. 

As we did for the no-vortex solution, we now intro- 
duce corrections from the edges. The improved solution 

{<t>Y ><f>? } is given by 



(41) 



where i = 1, . . . , ceil(N/2) and there is an additional 

4>N + l when N is even and 4>% + i = when N is odd. 
The double cell is used to simplify the calculations, as 
suggested by the spatial periodicity of the infinite-ladder 
solution. 

Again, far away from the ends the corrections are 
small and we linearize the governing equations (ft5l)-(ll6|) 
around the basic solution (|34|). Thus we obtain the fol- 
lowing system of coupled difference equations: 




Q -D l + 2E l = 
A - Ci+i +2F,=0 
-Ei cos (j)" d J= -Fi cos (f>" v J + A cos (j>X ve J 



(42) 
(43) 

(44) 

-Fi cos<pf ven = -Ei+i cos<^ d + C i+1 cos^ dd . (45) 

Eliminating Ei , F, and A we get a second-order differ- 
ence equation for C^: 



c i+2 + 2 7 a+i + q = o, 

with 7(7, /) given by 



(46) 



sin 7r/ — cb 2 



(47) 



(Here sa, sb, ca and cb are shorthand for sin 2na, sin 2irb, 
cos27ra and cos27r& respectively). This difference equa- 
tion has the general solution 



Ci = Pr 



■Qr2 



with r ff = -7 + y/j 2 - 1 ee e 1/Atf . 



(48) 
(49) 



And rff is related again to another characteristic pene- 
tration depth for the perturbations from the edges to die 
off, this time inside a region with the fully-frustrated so- 
lution. The coefficients P and Q have to be calculated 
numerically using the boundary conditions from nodes 1 
and N + 1. The spatial dependence of A is also of the 
same form Di = Rr l s + Sr^ 1 . Note that in both equa- 
tions, i is the number that indexes the double cell. 

Figure ^ (c) illustrates the accuracy of these approx- 
imate formulas. Specifically, we plot the spatial depen- 
dence of the predicted deviations 



Cf 



where 



D? 



is obtained from the numerical solution and 



i.vt 



is the infinite-ladder approximation. The expected 
exponential decay close to the edges with characteristic 
length Aff = 2/ In rff is verified in the figure. 



IV. DYNAMICAL DEPINNING TRANSITIONS 
OF THE ARRAY 

In this section, we describe the depinning transitions as 
seen in dynamical simulations of the ladder array at zero 
temperature. In the following section, we will explain 
these dynamical results by relating them to the bifurca- 
tions of the no- vortex, single- vortex, and fully-frustrated 
superconducting solutions. The main goal is to give a 
rigorous mathematical explanation of the following ob- 
servations: Dynamical simulations show that the array 
is superconducting at low values of the driving current 
/. As / is increased, the array remains superconducting 
until a critical current idcp(/) is reached, after which the 
array depins and develops a non-zero average voltage. 
This depinning current idep decreases monotonically as 
the frustration / increases from to 1/2. 

These numerical observations are all at the averaged 
level of the IV characteristics. They do not tell us any- 
thing about the detailed configuration of the individual 
junctions. In particular, there are several distinct super- 
conducting states (e.g., the nv, sv, and ff states discussed 
in Section III, and states containing multiple vortices) 
but these are indistinguishable on the IV curve. This 
ambiguity raises the question: what is the state of the 
ladder just before it depins? 
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We will show in Section [v| that for most values of /, the 
depinning of the ladder is caused by the destruction of the 
no- vortex state. Specifically, the global depinning current 
Idepif) can be predicted by calculating the current at 
which the no- vortex state is annihilated in a saddle-node 
bifurcation. The only exception occurs for values of / 
close to 1/2, where the depinning is due to saddle-node 
bifurcations of states of the fully-frustrated type. The 
noteworthy point here is that no other superconducting 
states play a role in the global depinning of the array. 
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FIG. 5. Snapshots of the time-evolution of the single-vortex 
solution for a 25 x 1 ladder array with f3 c = 10 and / = 0.2 

and increasing I. Initial (□), intermediate( ), and final 

(•) configurations are shown. At I — 0.134 < /lat the ini- 
tial configuration with a 27r-jump relaxes to the stationary 
single- vortex solution. When / = 0.136 ~ /lat the vortex 
becomes dislodged from the center. It moves slowly to the 
left and then stops at some intermediate position between the 
center and the edge. For Jlat < I = 0.2 < 7i c f t the vortex 
moves until it gets pinned near the edge, where the potential 
barrier is larger. Finally, at Ji c f t < I = 0.25 < Jdc P the vor- 
tex is expelled from the array and the no-vortex solution is 
recovered. 

However, the question arises as to how the depinning 
behavior would change when configurations with vortices 
are used as initial conditions (rather than the random or 
zero-phase initial conditions that we ordinarily use in our 



simulations of the IV characteristics). To address this 
issue, we perform dynamical simulations from an initial 
condition with one 27r-step in the middle cell of the array: 



(j)f(t = 0) = 0, 4>Y(t = 0) = 2tt 6 ( j - ceil 



N + 1 



where Q(x) is the Heaviside step function. This initial 
condition is not a solution of the system and, thus, under 
the dynamical equations (|l2|)-([l4|) it relaxes onto a true 
solution for the ladder. For most {/,/}, the single-vortex 
superconducting state (top panel of Figure |5|) is reached, 
i.e., a static configuration with a topological vortex in cell 
a of the array such that n a = 1 and rij = 0, Wj ^ a, j = 
1, . . . , N. For some ranges of / and /, this configuration 
is not dynamically stable and other solutions are found, 
as discussed in Section [v|. 

The numerical observations shown in Figure ^| depict 
the dynamical behavior of the single- vortex state for most 
values of /, i.e. approximately 0.12 < / < 0.37. They 
can be summarized as follows: For / > / m i n and small 
driving current /, the system relaxes onto a static single- 
vortex solution with a vortex in the middle of the array. 
The solution is slightly distorted as the current is in- 
creased, until at I = Jlat(/) the vortex moves from the 
center cell toward the left. (This current is analogous 
to the well-known Lobb-Abraham-Tinkham (LAT) de- 
pinning, current for two-dimensional Josephson junction 
arraysa). For currents very close to /lat, the vortex 
moves slowly and gets trapped in another cell somewhere 
between the center and the edge. For somewhat larger 
I > /lat, the vortex moves all the way to the left edge 
where it becomes pinned until, at a second critical cur- 
rent /left, it is expelled from the array and the no- vortex 
solution is recovered again. The no-vortex configuration 
then remains stable until, at idep, global depinning of the 
array occurs. If instead of placing the vortex in the mid- 
dle, we locate it closest to the right edge, it depins at a 
current /right < /lat, and moves toward the center. 

These observations can be clarified with the usual anal- 
ogjcj of the vortex as a damped particle moving in a 
sinusoidal potential under the action of a — x Lorentz- 
like force directly proportional to /. The maxima of the 
potential correspond to the vertical junctions and the 
minima are situated in the middle of the cells. Thus, 
the initial barrier which has to be overcome to begin the 
motion explains the critical /lat- Moreover, the open 
boundaries can be thought to produce an exponentially- 
decaying envelope superimposed on the sinusoidal poten- 
tial. Thus the motion of the vortex is easier when the 
vortex is close to the right edge and becomes increasingly 
difficult as the left boundary is approached. When /i e f t 
is reached, the vortex is able to overcome the edge bar- 
rier and is expelled from the array. Then the no-vortex 
configuration is recovered and no new vortex enters the 
ladder. 

We mentioned before that the described behavior is 
observed for values of the frustration contained between 







two limiting values. First, there is a minimum frustration 
/min ~ 0.12 below which the single- vortex solution never 
ensues from this initial condition; in fact, the system set- 
tles on the no-vortex superconducting solution. Second, 
for / larger than ~ 0.37, the vortex is not expelled from 
the array at the left edge before depinning. Dcpinning 
occurs in that case when vortices enter the array from 
the right edge. These observations will be clarified in 
Section [v|. 

We have also performed dynamical simulations for the 
multi-vortex case and reached similar conclusions for 
most values of /. In that case, the initial condition con- 
sists of Nf equally spaced 2-7r-steps which are then al- 
lowed to evolve dynamically. For I — and large enough 
/, the initial condition relaxes onto the expected solution 
with Nf vortices in the array. As the current is increased, 
these move towards the left end where they accumulate 
until they are expelled one by one at different currents. 
After this, the no-vortex solution is again recovered. 

On the other hand, the picture changes when / is close 
to 1/2. There, solutions of the fully-frustrated type are 
obtained from the multi- vortex initial condition and there 
is no expulsion of vortices from the ladder. Instead, the 
array depins globally at the current where the ff solution 
ceases to exist, i.e. Zdep(/ — ► 1/2) w Is,th- 

We conclude that in the ladder, the depinning of 
one vortex (or a train of vortices) is not equivalent to 
the global edge-dominated depinning of the device. As 
we will confirm in Section [v|, the no-vortex and fully- 
frustrated states are the relevant solutions for the de- 
pinning of the array; for moderate /, even if the initial 
conditions contain vortices, these are expelled from the 
array as the current is increased and, eventually, the no- 
vortex solution is recovered. For / close to 1/2, the sys- 
tem settles onto fully-frustrated solutions with distinct 
depinning properties. 

We have also identified three other critical currents 
related to the single- vortex configuration: /right) at which 
a vortex at the right edge begins to move; Jlat, at which 
dynamical depinning occurs for a single vortex centered 
in the middle of the ladder; and iieft, at which the vortex 
is expelled at the edge. 

All of these dynamical observations are explained in 
detail in the next section where they are compared to 
their exact mathematical descriptions. 



V. BIFURCATION ANALYSIS OF THE 
DEPINNING TRANSITIONS 

In this section, we use bifurcation theory to obtain 
exact criteria for all the critical currents of the no- 
vortex, single-vortex, and fully-frustrated states. We 
have checked consistently that these bifurcations explain 
the dynamical dep inn ing behavior of the array as de- 
Furthermore, analytical simpli- 



scribed in Section IV 



fications to some of those criteria will be deduced from 
approximations of the exact depinning results. 



The depinning of the ladder can be explained in dy- 
namical terms as follows: The linear stability of the su- 
perconducting states of the ladder as a function of / can 
be deduced from the Jacobian matrix Jd yn of the dynami- 
cal system (12)-(|l4|) for a given value of /, and for a given 
fixed point — in particular, the no- vortex, single- vortex or 
fully-frustrated state. If all the eigenvalues have negative 
real parts, the fixed point is linearly stable. As we in- 
crease /, some of the eigenvalues move to the right in the 
complex plane, and the fixed point becomes less stable. 
The critical current for a given fixed point is defined by 
the condition that the maximum of the real parts of the 
eigenvalues becomes positive. Then, to predict the global 
depinning current, we compare the critical currents of 
the different superconducting states, and take the max- 
imum of those. In other words, we predict that global 
depinning occurs when the "last" stable-state bifurcates. 

Recall that there are several scenariost2l by which a sta- 
ble fixed point can undergo such a bifurcation. First, in 
a zero-eigenvalue bifurcation, a single eigenvalue moves 
along the real axis, and passes from negative to posi- 
tive at the bifurcation. There are three main subtypes 
of zero-eigenvalue bifurcation: saddle-node, transcritical, 
and pitchfork. In the saddle-node bifurcation, a stable 
fixed point collides with a saddle point, and both are an- 
nihilated. In contrast, in the transcritical and pitchfork 
bifurcations, the stable fixed point is not destroyed - 
it continues to exist but goes unstable. A second sce- 
nario is provided by the Hopf bifurcation which involves 
a pair of complex conjugate eigenvalues passing through 
the imaginary axis from the left half plane to the right 
half plane - again this bifurcation destabilizes the fixed 
point, but does not destroy it. 

Since it can be showp-that Hopf bifurcations are not 
possible in this systemEU, we can simplify our calcula- 
tions by using the static system (fL5|)-(16) to identify the 
location of the zero-eigenvalue bifurcations. Those bifur- 
cation points are characterized by a change in the num- 
ber of fixed points and, thus, from the implicit function 
theoremES, the Jacobian matrix J of the static system 
has zero determinant there. Hence, we use the super- 
conducting (static) up-down symmetric system given by 
f(x)=0, 

fj(x) = I + sinxjv+i - sinxjv+i+i - sinxi, 

i=l,...,iV+l (50) 
fjv+i+i(x) = Xj - x j;+ i - 2xat +1+j + 27r(/ - rii), 

i=l,...,N (51) 

with x = (xi, . . . ,x 2i v+i) = {4>\ , ■ ■ ■ ,0w + i,0f > • ■ -i<I>n)- 
For a given / and a given superconducting state, we com- 
pute the bifurcating fixed point x* and its associated crit- 
ical current /*(/) at which det(J) = 0. To this end, we 
define an augmented algebraic system with the current / 
as an extra variable, and the constraint on the determi- 
nant as an extra equation. Then x* and I* are obtained 
by solving F(X*) = where X = (x, I) and 
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the depinning for the no-edge nv state occurs through a 
saddle-node bifurcation at J^ v th (/) = 1, V/. 
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FIG. 7. Saddle-node bifurcation of the no-vortex solution 
in a 15 x 1 ladder, (a) Value of the maximum of the real 
parts of the eigenvalues of the dynamic Jacobian matrix Jdyn 
for the stable (S, solid line) and unstable (U, dashed line) 
branches with / = 0.2 and (5 C = 0.25. They collide and 
annihilate at in a saddle-node bifurcation. We remark 
that, although the eigenvalues change with /3 C , the bifurcation 
point J* v is independent of /3 C . (b) Phase configuration of the 
vertical (left) and horizontal (right) junctions for the unstable 
branch (U) at / = 0.2 and 1 = 0. (c) Same as (b) for the 
stable branch (S). Note that both configurations (b)-(c) are 
up-down symmetric, and the unstable branch (b) corresponds 
to an inverted pendulum at the rightmost cell. 

To capture the observed /-dependence of the critical 
current I* v , we need a more careful approximation. We 
now present two such approximations to I* v that clarify 
the physical picture of the transition. 

The first strategy is to use the improved approxima- 
tion ( ]l8| ) for the no-vortex state (in which the edge ef- 
fects are included perturbatively), and then write down 
a simplified augmented system F^fX*) = for this so- 
lution, similar to the expression (p2])-(^3|) for the full 
(2N + 2)-dimensional system. We can then calculate the 
critical current I nv ,rcd for this reduced model, defined 
as the value of / where the perturbative solution under- 



goes a saddle-node bifurcation. This renders the calcu- 
lation analytically tractable since, for a given /, only 
three variables {P,Q, 1} suffice to describe the pertur- 
bative solution ( |l8|) — instead of I and 2N + 1 phases 
for the full solution. The simplified augmented system 
is constituted by (p8|), ( p9| ) together with the condition 
that the determinant of the 2x2 Jacobian matrix equals 
zero, corresponding to the zero-eigenvalue condition at a 
saddle-node bifurcation. 

Furthermore, since equations (|2^) and (|2^) arc un- 
coupled when N is not very small, even this three- 
dimensional system can be further reduced to a two- 
dimensional system with unknowns P and /: 



F \ = I + S m{nf-^(l-l 



~dP 



sin(arcsinI + P) (54) 



(55) 



where r = r(f,I) is given by (|25|). Note that this is a 
set of local equations referred to the rightmost end of the 
array. We numerically solve the 2x2 reduced system 
F^(P re d, Aiv,rod) — to obtain the approximate depin- 
ning current / n v,red(/) and the value of the rightmost 
phase P rc d at the bifurcation. 

Figure || shows that / nv , ro d(/) predicts the exact J£ v (/) 
reasonably well. As expected, the prediction gets worse 
as / nears 1/2 since the perturbative approximation of 
the no- vortex solution is less accurate in that limit. 

Both the analyses of the eigenfunctions of the full 
2iV + 2 system (|g)-(|3|) and of the reduced system (||)- 
( [55| ) indicate that the global depinning of the ladder is 
caused by a local instability of the rightmost junction 
of the array. This is consistent with physical arguments 
which imply that after depinning occurs, vortices propa- 
gate in the array in the — x direction under the effect of 
a magnetic Magnus (Lorentz-like) force. 

The key role played by the rightmost junction suggests 
a second simplification, which we call a heuristic crite- 
rion for depinning. This criterion connects the global 
depinning of the ladder with the much simpler depin- 
ning transition in a single junction. Recall that when 
the phase of a single junction reaches 7r/2, its supercon- 
ducting solution is destroyed in a saddle-node bifurca- 
tionE]. Therefore, we intuitively propose that when the 
phase of the rightmost junction reaches 7r/2, the ladder 
depins. Replacing the no-vortex solution by its pertur- 
bative approximation (|l8|), we solve for the current I n v,th 
by imposing 



%+i = t/2 



which implies 



arcsin/ nVj th + P(I n 



Ah) 



ir/2. 



(56) 



(57) 



Then, from (p9|), we obtain an implicit transcendental 
equation for J nv ,th(/): 
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FIG. 10. Saddle-node bifurcation of the single-vortex solu- 
tion at the center of a 15 x 1 ladder, (a) Value of the maximum 
of the real parts of the eigenvalues of Jd yn for the stable (S, 
solid line) and unstable (U, dashed line) branches with / = 0.2 
and (i c = 0.25. Jg V is the point where a saddle-node bifurca- 
tion occurs for this particular configuration. As in Fig. (t| 
the bifurcation point is independent of /3 C . (b) Phases of the 
vertical (left) and horizontal (right) junctions for the unsta- 
ble branch (U) at / = 0.2 and 1 = 0. (c) Same as (b) for 
the stable branch (S). Again, both configurations (b)-(c) are 
up-down symmetric. Note that the stable branch (c) is asso- 
ciated with a vortex located at the center of a plaquette (a 
minimum of the potential), whereas the unstable branch (b) 
corresponds to a vortex centered around a junction (a local 
maximum of the potential energy). 



2. Symmetry-breaking bifurcation of the sv solution 

We noted above the numerical observation!!!! that, 
when performing dynamical simulations, the static 
single-vortex solution is unstable below a critical field 
/min(-0- We show now that this is the result of a 
symmetry-breaking instability which is mathematically 
related to another zero-eigenvalue bifurcation. Therefore, 
once again, the dynamical / m m coincides with a static f* v 
calculated from the augmented system (p2|)-(|53|) as the 
value of / where the determinant of the static Jacobian 
matrix J is zero and, thus, a change in the number of 



fixed points is likely. Excellent agreement between / m ; n 
and f* v is shown in Figure ^. 

Figure [ll] depicts detailed information about this bi- 
furcation. Specifically, it shows three single- vortex states 
that co-exist for / slightly greater than f* v . These states 
appear very similar, but on close inspection, one notices 
that two of the states are asymmetric: cj)^ ^ — </>^ (this 
is especially clear for the central plaquette j — a.). As 
/ ~^ /sv from above, these asymmetric states - which 
are unstable - simultaneously collide with the stable sym- 
metric state, rendering it unstable. To visualize this tran- 
sition in greater detail, Figure [ll] plots the asymmetry for 
the central plaquette, S = (j)^ + (f^ , as a function of the 
frustration /. The symmetric state exists both above and 
below the bifurcation, and satisfies S — throughout. 
The two unstable branches join the symmetric branch at 
/ = fsv ■ The scenario depicted in Figure [ll], common in 
symmetric jsv stems, is known as a subcritical pitchfork 
bifurcatiorJiS. 

From this, we conclude that there exists a region of the 
(/, I) plane (Figure ^) where the single- vortex configura- 
tion exists but is always dynamically unstable. In this re- 
gion, the vortex (magnetic flux) is expelled transversally 
from the array through transient modes which do not pre- 
serve the up-down symmetry of the horizontal phases. 

We can also derive an analytical expression / sv ,th for 
the critical field f* v . As given in (11), all up-down sym- 
metric fixed points become unstable when the phase of 
any of the horizontal junctions is larger than 7r/2 in abso- 
lute value. Since this largest phase occurs at the central 
plaquette (as seen in Figure [ll]) , the following stability 
criterion ensues from the condition <f>^ = —ir/2: 



27r/ SVjt h + [P(/sv,th) — Q (/sv,th)] = 7T, 



(59) 



where we have used (|3(j) and (Q) and we consider a long 
array such that the effect of the edges on the middle cell 
can be neglected. Here P, Q' have to be calculated from 
the boundary conditions (current conservation) at nodes 
a and a + 1 : 

J + sin(7r/-P(r-l)/(2r)) = 
sin(arcsin i+P) - sin(7r/ - (Q' - P)/2) 

I - sin(7r/ - (Q' - P)/2) = 
sin(arcsin/+Q') + sin(7r/ + Q'(r - l)/(2r)) (60) 



with r given by (25). Numerical solution of the critical 
condition ( |59| ) yields the curve / sv ,th(-0 which is almost 
indistinguishable from the curve for f* v from the full aug- 
mented system (Figure |9|). 

It is especially interesting to check the case 1 = 0. 
Then P = —Q' and the critical condition ([59|)-(|60"|) sim- 
plifies to 



sv,th 



f" 



sv.th 



7r(r — 1) 
2~r 



= cos(7r/ s ° V)th ) - 1, 
(61) 
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which can be solved numerically to give f° th = 
fsv,th(I = 0) = 0.1193, in agreement with numerical find- 
ings from previous dynamical simulationsEj. 
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FIG. 11. Subcritical pitchfork bifurcation of the sin- 
gle-vortex solution at / a * v . (a) Value of the up-down asym- 
metry S = 4>a + 0a for the two asymmetric (AS + , AS - ) and 
symmetric (S°) branches with 1 = 0. The calculation of the 
dynamical eigenvalues (not shown) indicates that the stable 
symmetric branch becomes unstable at / s * v when it collides 
with the two unstable asymmetric branches. For / < /* v , only 
the unstable symmetric branch survives, (b) Phases of the 
vertical (left) and horizontal (right) junctions for the asym- 
metric branch (AS + ) at / = 0.14 and 7 = 0. (c) and (d) are 
the same as (b) but for the S° and the AS - branches respec- 
tively. Note that although all three vertical configurations on 
the left panels look very similar, the (b) and (d) AS horizontal 
configurations on the right are up-down asymmetric. 



3. Single-vortex configurations at the 



So far we have concentrated on a particular single- 
vortex state, namely one where the vortex occurs in the 
middle cell of the array. But there are many other single- 
vortex states, each differing from the previous one by dis- 
placing the vortex by one cell to the right or to the left. 
Each of those configurations becomes unstable through 
bifurcations similar to those discussed above for the case 



of the centered single-vortex state. For most values of 
/, as discussed in Section IV, when the driving current 
I is increased, a vortex in the ladder moves to the left 
(getting pinned in cells closer to the left edge as I grows) 
until it is expelled from the array at the left boundary. 
Thus, in explaining the effect of the edges on the sv solu- 
tion we are most interested in two critical currents: the 
critical current / r i g ht , at which a vortex at the right edge 
begins to move, and the critical current ii e f t , at which 
the vortex is expelled at the left edge of the ladder. 

To predict / r i g ht , we analyze the sv configuration with 
the vortex placed in the rightmost cell (the 000 ... 01 
configuration). The results of the analysis are similar 
to those for the vortex in the center. As shown in Fig- 
ure |l^(a), this state can cease to exist through a saddle- 
node bifurcation (J* v right ), or become unstable through a 



symmetry-breaking pitchfork bifurcation (/* v 



The 



agreement with the dynamical simulations is excellent. 

The rigorous explanation of Ji e £t turns out to be slightly 
more complicated. A careful examination of the numerics 
reveals that, depending on the value of /, the vortex can 
be expelled from the array in one of two ways: from the 
leftmost cell (cell number 1) as expected, or directly from 
cell number 2. Thus, we need to examine the dynamical 
stability of two sv states: the one with the vortex in cell 
1 (1000 ... 00) and in cell 2 (0100 . . . 00). 

The dynamic stability of the 0100 ... 00 configuration 
contains no new elements. The two observed bifurca- 
tions (a saddle-node I* v loft2 and a subcritical pitchfork 
fsv ieft2) are similar to those explained above. The re- 
sults are presented in Figure [l2](b) where the saddle- 
node bifurcation is seen to explain the dynamical 7i e f t 
for / < 0.29. 

However, it is the stability analysis of the 1000 ... 00 
state that explains the expulsion of the vortex at the edge 
for / > 0.29. As shown in Figure fl^(c), this configura- 
tion presents the usual saddle- node lj* v i eftl ) bifurcation. 
The pitchfork (f* v lcftl ) bifurcation is barely visible in the 
figure. There is also another saddle-node bifurcation at 
low / and high / which is irrelevant for the depinning 
considered here. 

The results of this section are summarized in Fig- 
ure |l2|(d), which indicates the region where at least one 
single- vortex configuration in a ladder array with edges is 
dynamically stable. Figur e |l2| (d) is, in essence, the union 
of Figure [| with Figure |12|(a)-(c) and shows how the 
sv solutions either cease to exist through a saddle-node 
bifurcation when / is increased for most values of the 
frustration /, or become unstable through a symmetry- 
breaking pitchfork bifurcation for small values of /. 

Within this picture, the dynamical behavior in Fig- 
ure^, where the vortex propagates along the ladder in the 
interval J* V)right (/) <I < I* v , lclt (f), is the result of a suc- 
cession of saddle-node bifurcations of single- vortex states 
situated in contiguous cells until the vortex is expelled 
at the left edge. This is in contrast with the symmetry- 
breaking exit of the vortex for / < f* v (I) where the flux 
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is expelled in the transversal direction. 



(a) 



v L. 




0.1 0.2 0.3 0.4 
frustration, f 



FIG. 12. (a)-(c) Stability diagrams for the sv configuration 
when the vortex is placed close to either edge. The symbols 
denote critical currents measured through dynamical simula- 
tions. The solid lines (l£ v ) correspond to saddle-node bifurca- 
tions and the dashed lines (fs V ) to subcritical pitchfork bifur- 
cations where the flux is expelled transversally. For instance, 
from (a) the vortex in the rightmost cell moves to the left 
for / > J* Vj right . And ii e ft corresponds to the expulsion of the 
vortex at the left edge for / < 0.37. In the region delimited by 
bright and Zioft, there is a cascade of saddle-node bifurcations 
which leads to the behavior observed in Figure [s| For the top 
three panels: (a) corresponds to the vortex in the rightmost 
cell; (b) in the second cell; and (c) in the first (leftmost) cell, 
(d) shows the region where at least one sv configuration (any- 
where in the ladder) is stable. It is the union of panels (a)-(c) 
and Figure ^| 

This description is valid for / < 0.37. However, a 
subtle variation is observed beyond that point for the sv 
solutions. The comparison of Figure |l2|(c) and Figure || 
shows that the depinning of the sv and nv states is simi- 
lar for high /, i.e. I* v loft « I* v , for / > 0.37. This is due 
to the fact that at high / the vortex is not expelled from 



the array before depinning. Instead, the sv configuration 
with the vortex at the left edge ceases to exist through 
a saddle-node bifurcation which, similarly to what hap- 
pens to the nv state, is localized on the right edge of 
the ladder. The depinning transition of the single-vortex 
state at high f can be described in a similar fashion as 
the depinning of the no-vortex configuration: by the nu- 
cleation of vortices on the right edge of the ladder. Thus, 
the single- vortex configuration plays no role in the global 
depinning of the array: only the behavior of the nv and 
ff states has to be considered. 



C. Bifurcations of the fully- frustrated state 

The depinning of the ladder array as / — > 1/2 is deter- 
mined by the stability of the fully-frustrated configura- 
tion. This can be readily seen by inspection of Figure ^ 
which shows that at high / the depinning of the no- vortex 
state occurs before the global depinning of the array. 

To clarify the importance of these transitions we follow 
the same scheme as above. We first calculate the zero- 
eigenvalue bifurcations of this state and show that the 
depinning of the ladder at / ~ 1/2 is indeed explained 
by a saddle-node bifurcation of the ff state. Then we 
obtain analytical approximations ("th") to the exact bi- 
furcations, using stability criteria for the infinite-ladder 
ff configuration. 

In principle, the characterization of the ff bifurcations 
is more intricate than for the nv and sv configurations 
above since in the finite ladder several states of the ff 
type could play a role in the depinning. First, there ex- 
ist different states in ladders with odd and even number 
of cells, as seen when the . . . 010101 . . . alternating vor- 
tex pattern is fitted into a finite length. Second, there 
are many ff-states very close energetically with different 
dynamic stability structure. However, we will show that 
the landscape of relevant solutions is indeed clear, and 
depinning can be assigned to instabilities of one of those 
configurations. 

Consider first an even ladder with ten plaquettes as an 
example. Of the several states of the ff type, three can be 
thought as relevant: (la) 1010101010, (lb) 0101010101, 
and (lc) 1101010100. We have analyzed the stability of 
these three states and conclude that only configuration 
(la) is relevant for depinning of the even ladder. Both 
(lb) and (lc) cease to exist or become unstable at lower 
critical currents. This is physically reasonable since (lb) 
will tend to move one cell to the left under the action of a 
driving current to produce (la), while (lc) ceases to exist 
through a low current saddle-node bifurcation caused by 
the expulsion of the vortex at the left boundary. 

The case of the odd ladder with N plaquettes has one 
further complication, namely that we cannot have ex- 
actly N/2 vortices in the array. Thus, the "pure" ff state 
(as calculated for the infinite case) is not possible un- 
der these topological constraints. However, states simi- 
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Our solution resembles the sine-Gordon kink in that it 
describes a 27r-jump with odd symmetry with respect to 
xq = a + 1/2. However, both the functional forms and 
the characteristic lengths are different. Figure [n] shows 
that the nume rics are better approximated over a wi de 
range of / by (A2) than by the sine-Gordon kink (|Al|) . 

Given the relative inappropriateness of the sine- 
Gordon kink as a model of the vortex in the no- 
inductance ladder array, we conclude that even for the 
static case it is an oversimplification to reduce the ladder 
to a one-dimensional parallel array where the horizontal 
junctions are approximated by an effective inductance. 
Other dynamic phenomena observed in the ladder, e.g, 
the dynamical mechanism of retrapping from the whirling 
mode^-seinforcc this statement and will be discussed else- 
whereEil. 



APPENDIX B: NO- VORTEX SOLUTION AND 
DEPINNING WITH SELF-INDUCTANCE 

Following on the concepts and notation in Sec- 



tions III and IV , we briefly consider the ladder with self- 
inductance. This is a first approximation to explore the 
effect of self-fields on the no-vortex solution and, conse- 
quently, on the global dcpinning current of the array. 

When self-inductances are included, the time- 
independent governing equations (|l5|)-(|l6|) become 



I + sm(j)f_ 1 



sin (f>? 



20? = 2 7 r(n j - /) - L"VA 2 , 



(Bl) 
(B2) 



The procedure to account for the open boundaries is 
identical to that presented in Section III . The results are 
also similar: the corrections decay exponentially from the 
edges as in (|3|)-(|2^). However, the characteristic length 
1/ In r s is now given by r s = a s + \J a 2 s — 1 where 



A 



a. = 1 



cos 



LH* 



1 



1 

2A 2 



cos 



lH* 



(B5) 
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FIG. 17. Effect of the self-inductance on the depinning cur- 
rent. We show how the heuristic depinning approximation 
-fnv,th(/), obtained from (B6), varies with the self-inductance 
A = Lj I L a - Note that the rest of the paper deals with the 
limiting case A = oo where all inductances are neglected. 



where i™ is the mesh current in plaquette j and A 2 = 
Lj/L s is the two-dimensional penetration depth and a 
measure of the discreteness of the array. Due to the geo- 
metrical constraints of the ladder, it is readily seen that 
T" 1 = —I? in this case. 

For the no-vortex solution we still have (j)J 
and *h „ j.h _ j,h* 



far from the 



edges. 



From (Bl) we then get 



arcsin/, 



(B3) 



as for the case with no inductance (A = oo). However, 
the horizontal phase is different and has to be calculated 
from the following nonlinear equation: 



1 

2A 2 



• J.H* /• 



(B4) 



This is a particular-ease of Kepler's equation, studied in 
celestial mechanicsEll, which can be solved through the 
method of successive approximations. In our case, a good 
approximation to <f> H over the whole range of / is given 
by the first iteration of that method as: 



, 1 . f 2A 2 

7T f TTT Sin -tt 

J 2A 2 ll + 2A 2 



To assess how the inductance affects the depinning cur- 
rent, we use (B5) and the <fi H obtained from (B4) with 
a heuristic criterion similar to (58), 



arcsin(l - 7 nv ,th) 



• H* 



2r, 



■ arccos L 



nv,th 7 



(B6) 



to calculate / n v,th(/, A 2 ). The results in Figure [l7] show 
that, for fixed A, 7 n v,th(/) still decreases monotonically 
with /. However, for a given /, the value of the depin- 
ning current increases as A diminishes. This is expected 
on physical grounds since an increase in the inductance 
L s implies larger self-fields which oppose the external ap- 
plied magnetic field, thus decreasing the effective mag- 
netic flux through the array. 
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Erratum: Superconducting states and 
depinning transitions of Josephson lad- 
ders [Phys. Rev. B 57, 1181 (1998)] 

Mauricio Barahona, Steven H. Strogatz, and 
Terry P. Orlando 

The first formula of Section VI should include 
a contribution from the driving current. In our 
choice of gauge, 

N+l 

E = -Y,<xx<l>j- I Y, fl- 
an j=l 

Accordingly, Eq. (64) becomes 

Vl - P+2 cos(7t/ NV -ff)+/ arcsin(7)-2 sin(7r/ NV -FF) (1-2I 2 / L)~ 1/2 -I arcsin( v /Z/2) = 0, 

and /nv-ff (I = 0) = 1/3 remains unaltered. 
Figure 15 is very slightly modified as shown. The 
conclusions do not change. 




FIG. 18. (Corrected Fig. 15) Phase diagram of the NV, 
SV and FF solutions in the infinite (no edges) ladder — also 
relevant for the ring ladder. Note that region H shrinks as the 
number of junctions is increased (i.e. / sv -ff tends asymptoti- 
cally to /nv-ff as N — » oo). Thus, the region H in the figure 
is only indicative. 



23 



